BuildReachNetwork Subroutine

SUBROUTINE BuildReachNetwork &
!
(maxReachLength, slopeCorrection, flowDirection, flowAcc, dem, &
    fileExport, shpExport, streamNetwork )


IMPLICIT NONE

!arguments with intent(in):
REAL (KIND = float), INTENT (IN) :: maxReachLength  !!max length of a reach (m)
REAL (KIND = float), INTENT (IN) :: slopeCorrection !! slope value to correct negative values
TYPE (grid_integer), INTENT (IN) :: flowDirection !! flow direction
TYPE (grid_integer), INTENT (IN) :: flowAcc !! flow accumulation
TYPE (grid_real),    INTENT (IN) :: dem !! digital elevation model
INTEGER (KIND = short), INTENT (IN) :: fileExport, shpExport

!arguments with intent (out):
TYPE (ReachNetwork), INTENT (OUT) :: streamNetwork

!local declarations
LOGICAL                         :: outletSection
LOGICAL                         :: endOfReach
INTEGER                         :: countReaches
INTEGER                         :: maxOrder
INTEGER                         :: i,j,l !loop index
INTEGER                         :: row,col !current cell
INTEGER                         :: iDown, jDown !downstream cell
INTEGER                         :: n_cells
REAL                            :: reachLength

TYPE (grid_integer)             :: orderBeginning
TYPE (grid_integer)             :: orders
TYPE (grid_integer)             :: split
TYPE (Coordinate)               :: point1, point2

!- -----------------------------end of declarations----------------------------

!------------------------------------------------------------------------------
!(1.0) Initialize orders, orderBeginning and split grids
!------------------------------------------------------------------------------

CALL NewGrid (orders, flowDirection)

CALL NewGrid (orderBeginning, flowDirection)

CALL NewGrid (split, flowDirection)

!------------------------------------------------------------------------------
!(2.0) Build Horton orders grid
!------------------------------------------------------------------------------

CALL HortonOrders (flowDirection, orders, maxOrder)

!------------------------------------------------------------------------------
!(3.0) Build orderBeginning and split network
!------------------------------------------------------------------------------

CALL MarkOrderBeginning (orders, flowDirection, orderBeginning)

!split network
CALL SplitNetwork (orders, flowDirection, split)

!------------------------------------------------------------------------------
!(4.0) Build reaches
!------------------------------------------------------------------------------

countReaches = 0
NULLIFY(list)  ! at the beginning list is empty

!initialize points to compute reach length
point1 % system = flowDirection % grid_mapping
point2 % system = flowDirection % grid_mapping

DO l =1, maxOrder

    CALL Catch ('info', 'RiverDrainage', 'Elaborating reaches of stream order: ', &
                 argument = ToString(l))


    DO j=1,orderBeginning%jdim         ! scan orderBeginning grid
      DO i=1,orderBeginning%idim

        IF(orderBeginning%mat(i,j) == l) THEN
           countReaches    = countReaches + 1 
           reachLength     = 0.
           n_cells         = 0
           row             = i
           col             = j
           endOfReach      = .FALSE.
           outletSection   = .FALSE.

           IF(.NOT.ASSOCIATED(list)) THEN
               ALLOCATE(list)     !set first reach
               current => list
           ELSE
               ALLOCATE(current%next) !allocate next reach
               NULLIFY(current % next % next)
               current => current%next
           ENDIF

           current % i0    = i
           current % j0    = j
           current % order = l
           current % id    = countReaches

    	  
           !follow the reach until an order change or to the outlet, 
           !if it is the maximum order    
           DO WHILE (.NOT. endOfReach) 
           
          
             CALL DownstreamCell (row,col,flowDirection%mat(row,col), &
                                 iDown,jDown)
             n_cells      = n_cells + 1
             
             CALL GetXY (row, col, flowDirection, point1 % easting, point1 % northing)
             CALL GetXY (iDown, jDown, flowDirection, point2 % easting, point2 % northing)
             reachLength  = reachLength + Distance (point1, point2) !length (m)
             
		     IF ( CheckOutlet (row,col,iDown,jDown,flowDirection) ) THEN
			    outletSection = .TRUE.
			    endOfReach    = .TRUE.
		     ENDIF

             IF ( outletSection ) THEN
               current % i1 =  row
               current % j1 =  col
		       current % n_cells = n_cells ! unit: cells (integer)
		       current % length = reachLength ! !length (m)
		       current % area = flowAcc % mat(row,col) ! * &
		                        !flowAccumulation % cellsize**2
             
		       		  		       		       		       		      	        		          
		       ! if the reach length exceedes the maximum length =>
		       ! break the reach
             ELSEIF (maxReachLength > 0. .AND. &  !if maximum length = 0 I do not break the reach
		             reachLength  >= maxReachLength ) THEN
		       current % i1 = row
		       current % j1 = col
		       current % n_cells = n_cells
		       current % length = reachLength
		       !current % routingMethod = default
		       current % area = flowAcc % mat(row,col) ! * &
		                        !flowAccumulation % cellsize**2 !(m2)
		       !set elements for the following reach
		       ALLOCATE(current % next)
		       NULLIFY(current % next % next)
		       current => current % next
		       current % i0    = row
		       current % j0    = col
		       current % order = l
		       countReaches  = countReaches + 1
		       reachLength   = 0. 
		       n_cells       = 0
		       current % id  = countReaches
		       
		     END IF	
		     
		     IF ( .NOT. IsOutOfGrid(iDown,jDown,orders)) THEN
		         IF (orders%mat(iDown,jDown) > l) THEN
                   endOfReach = .TRUE.
                   current % i1 = iDown
                   current % j1 = jDown
                   IF (n_cells == 0) THEN !this case occurs when a reach has just been terminated because too long
                       current % n_cells = 1
                       CALL GetXY (row, col, flowDirection, point1 % easting, point1 % northing)
                       CALL GetXY (iDown, jDown, flowDirection, point2 % easting, point2 % northing)
		               current % length = Distance (point1, point2)
                   ELSE
		               current % n_cells = n_cells
		               current % length = reachLength 
		           END IF
		           current % area = flowAcc % mat(row,col) ! * &
		                            !flowAccumulation % cellsize**2
		         END IF  
		     END IF
		      		 			                  
             row = iDown
             col = jDown
                         
           END DO                     
        
        ENDIF
                
      ENDDO
    ENDDO

ENDDO


!set other parameters
!scan all list elements
current => list
DO 

  !calculate slope from digital elevation model and check negative slope
  current % slope = ( dem % mat(current%i0,current%j0) - &
                      dem % mat(current%i1,current%j1) ) / &
                      current % length
                      
  IF ( current % slope <= 0. ) THEN
    IF ( slopeCorrection > 0. ) current % slope  = slopeCorrection
  ENDIF
    
   !calculate spatial coordinate
   CALL GetXY (current % i0, current % j0, dem, current % x0, current % y0)
   CALL GetXY (current % i1, current % j1, dem, current % x1, current % y1)
  
  IF ( .NOT. ASSOCIATED (current % next) )THEN !found last element of list
    EXIT
  END IF
  
  current => current % next
  
END DO


!------------------------------------------------------------------------------
!(5.0) remove temporary variables
!------------------------------------------------------------------------------

CALL GridDestroy ( orders )

CALL GridDestroy ( orderBeginning )


!------------------------------------------------------------------------------
!(6.0) export files
!------------------------------------------------------------------------------

IF ( fileExport > 0 ) THEN
    CALL ExportReaches ( )
END IF

IF ( shpExport > 0 ) THEN
    CALL ExportShape ( dem, flowDirection )
END IF

!------------------------------------------------------------------------------
!(7.0) populate stream network 
!------------------------------------------------------------------------------
ALLOCATE ( streamNetwork % branch (countReaches) ) 

streamNetwork % nreach = countReaches

current => list
DO i = 1, countReaches
    
  streamNetwork % branch (i) % id = current % id
  streamNetwork % branch (i) % i0 = current % i0
  streamNetwork % branch (i) % j0 = current % j0
  streamNetwork % branch (i) % i1 = current % i1
  streamNetwork % branch (i) % j1 = current % j1
  
  streamNetwork % branch (i) % x0 = current % x0
  streamNetwork % branch (i) % y0 = current % y0
  streamNetwork % branch (i) % x1 = current % x1
  streamNetwork % branch (i) % y1 = current % y1
  
  streamNetwork % branch (i) % ncells = current % n_cells
  streamNetwork % branch (i) % slope = current % slope
  streamNetwork % branch (i) % length = current % length
  streamNetwork % branch (i) % area = current % area
  streamNetwork % branch (i) % order = current % order
    
  current => current % next

END DO
  
RETURN
END SUBROUTINE BuildReachNetwork